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Abstract 



Dynamics of FitzHugh-Nagumo (FN) neuron ensembles with time-delayed 
couplings subject to white noises, has been studied by using both direct sim- 
ulations and a semi-analytical augmented moment method (AMM) which 
has been proposed in a preceding paper [H. Hasegawa, Phys. Rev E xx, 
yyyy (2004)]. For iV-unit FN neuron ensembles, AMM transforms original 
2A r -dimensional stochastic delay differential equations (SDDEs) to infinite- 
dimensional deterministic DEs for means and correlation functions of local 
and global variables. Infinite-order recursive DEs are terminated at the fi- 
nite level m in the level-m AMM (AMMra), yielding 8(m + l)-dimensional 
deterministic DEs. When a single spike is applied, the oscillation may be 
induced when parameters of coupling strength, delay, noise intensity and/or 
ensemble size are appropriate. Effects of these parameters on the emergence 
of the oscillation and on the synchronization in FN neuron ensembles have 
been studied. The synchronization shows the fluctuation-induced enhance- 
ment at the transition between non-oscillating and oscillating states. Results 
calculated by AMM5 are in fairly good agreement with those obtained by 
direct simulations. 
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I. INTRODUCTION 



There have been many studies on effects of noises in dynamical systems with delays. 
Complex behavior due to noise and delay is found in many systems such as biological systems, 
signal transmissions, electrical circuits and lasers. Systems with both noises and delay are 
commonly described by stochastic delay differential equations (SDDEs). In recent years, 
linear SDDEs of Langevin equation are beginning to gain much attention [1]- [6]. The 
parameter range for the stationary solutions of the Langevin equation has been examined 
with the use of the step by step method [1], the moment mothod [2] and the Fokker-Planck 
equation (FPE) method [3] [4]. 

When we pay our attention to living brains, various kinds of noises are reported to be 
ubiquitous. A study on noise effects has been one of major recent topics in neuronal systems. 
It has been shown that the response of neurons may be improved by background noises. The 
typical example is the stochastic resonance in which weak noises enhance the transmission of 
signals with the subthreshold level. The transmission delay is inherent because the speed of 
spikes propagating through axons is finite. Conduction velocity ranges from 20 to 60 m/s, 
leading to non-negligible transmission times from milliseconds to hundreds milliseconds. 
Although an importance of effects of delay has been not so recognized as that of noises, 
there is an increasing interest in the complex behavior of time delays, whose effects have 
been investigated by using integrate-and-fire (IF) [8]- [12], FitzHugh-Nagumo (FN) [13]- [15], 
Hindmarsh-Rose (HR) [16], and Hodgkin-Huxley (HH) models [10] [11] [17] [18]. Exposed 
behaviors due to time delays are the multistability and bifurcation leading to chaos. 

There are two difficulties in studying combined effects of noise and delay in brains. One is 
that the system is usually described by nonlinear SDDEs, which are generally more difficult 
than linear SDDEs. Dynamics of individual neurons includes a variety of voltage dependent 
ionic channels which can be described by nonlinear DEs of Hodgkin-Huxley-type models, 
or of reduced neuron models such as IF, FN and HR models. The other difficulty is that 
a small cluster of cortex consists of thousands of similar neurons. For a study of dynamics 
of noisy neuron ensembles with time-delayed couplings, we have to solve high-dimensional 
nonlinear SDDEs, which have been studied by direct simulations (DSs) [19] [20] and by 
analytical methods like FPE [21]. Simulations for large-scale neuron ensembles have been 
made mostly by using IF, FN, HR and phase models. Since the time to simulate networks 
by conventional methods grows as N 2 with N, the size of the ensemble, it is rather difficult 
to simulate realistic neuron clusters. Although FPE is a powerful method in dealing with 
the stochastic DE, a simple FPE application to SDDE fails because of its non-Markovian 
property [3] [5]. 

In a preceding paper [22] (which is referred to as I hereafter), the present author has 
developed an augmented moment method (AMM) for SDDE, employing a semi-analytical 
dynamical mean-field approximation (DMA) theory [23] [24]. In I, AMM is applied to an 
ensemble described by the delay Langevin model, transforming the original iV-dimensional 
SDDEs to infinite-dimensional DEs which are terminated at finite level m in the level- 
m AMM (AMMm). Model calculations in I with changing the level m have shown that 
calculated results converge at a fairly small m. Actually results obtained by AMM6 are 
in good agreement with those by DSs for linear and nonlinear Langevin ensembles. It 
has been demonstrated in I that AMM may be a useful tool in discussing dynamics and 
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synchronization of ensembles described by SDDEs. 

It is the purpose of the present paper to apply AMM to FN neuron ensembles with time- 
delayed couplings. In the next Sec. II, we apply our AMM theory to nonlinear SDDEs of 
Af-unit FN neuron ensembles, in order to get the infinite-dimensional deterministic DEs for 
the correlation functions of local and global variables. Infinite-dimensional recursive DEs 
are terminated at the finite level m in AMMm. In Sec. Ill we report model calculations, 
showing that results of our AMM are in good agreement with those of DSs. Section IV is 
devoted to conclusions and discussions. 

II. FN NEURON ENSEMBLE 

A. Adopted model and method 

Dynamics of a neuron ensemble consisting of A^-unit FN neurons (N > 2), is described 
by the 2A^-dimensional nonlinear SDDEs given by 

F[x u (t)] - cx 2i {t) + (-J-) £ w i3 G( Xlj (t - r y )) + at) + I {e) (t), (1) 

bx u {t) - dx 2i {t) + e, (i — l — N) (2) 

where F[x(t)\ = k x(t) [x(t) -h][l- x(t)\, k = 0.5, h = 0.1, b = 0.015, c = 1.0, d = 0.003 
and e = [23] [25], and xu and x 2i denote the fast (voltage) and slow (recovery) variables, 
respectively. The third term in Eq. (1) stands for interactions with the uniform couplings 
of Wij = w and delay times of r^- = r, and the sigmond function G(x) given by G(x) = 
l/(l+exp[—(x—9)/a\), 9 and a denoting the threshold and the width, respectively [26]. The 
all-to-all couplings have been widely employed in theoretical studies. The assumed constant 
delay may be justified in certain neural networks [27]. The fourth term of Eq. (1), £i(t), 
denotes the Gaussian white noise given by < >= and < &(t) >= f3 2 5 i:j 5(t - t') 
where (3 denotes the magnitudes of independent noises and the bracket < • > the stochastic 
average [28]. The last term in Eq. (1), I^ e \t), denotes an external input whose explicit form 
will be shown later [Eq. (31)]. 

We apply our AMM developed in I to FN neuron ensemble given by Eqs. (1) and (2), 
defining global variables for the ensemble given by 

= ^ E M*), « = 1, 2, (3) 

i 

and their averages by 

fM K {t) = < X K (t) > . (4) 
We define the correlation functions between local variables, given by 

7k,a(*> = ^7 E < 8x Ki (t) 5x Xi (t') >, k, X = 1, 2 (5) 



dx u (t) 
dt 

dx 2 j{t) 
dt 
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where Sx K i(t) = x K i(t) — p K (t). Similarly we define the correlation function between global 
variables, given by 



p K ,x(t,lf) = < 8X K (t) 5X x (t>) >, 



(6) 
(7) 



where SX K (t) = X K {t) — p K (t). Conventional variances and covariances are given by Eqs. 
(5)-(7) with t = t', for which the symmetry relations: 7i, 2 (M) — 7 2 ,i(M) and Pi >2 (t,t) = 
P2,i(t,t), are hold. It is noted that 7 /C)I/ (i, t) (k, A = 1, 2) expresses the spatial average of 
fluctuations in local variables of x Ki while p KtV (t,t) denotes fluctuations in global variables 
oiX K . 

After our previous studies [22-24], we have assumed that the noise intensity f3 is weak 
and that the distribution of state variables takes the Gaussian form concentrated near the 
means of (pi, p 2 ). The second assumption is justified from numerical calculations for single 
FN [29,30] and HH neurons [31,32]. We will obtain infinite-order equations of motions for 
means, variance and covariances defined by Eqs. (5)- (7). They will be terminated at the 
level m in AMMm. Readers who are not interested in mathematical details, may skip to 
Sec. IIC. 



B. Equations of motions 

After some manipulations, we get DEs for p K (t), 7 K ,„(i, t) and p KjV (t, t) (k, v — 1, 2) given 
by (for details see appendix A) 



dpi(t 



dt 
dp 2 (t 



dt 

d7i,i(M 



dt 

<*72,2(M 



dt 

dji,2(t,t 



dt 

dpi,i(t,t 



dt 

dp 2 ,2(t,t 



dt 

dpi, 2 (t,t 



dt 



with 



fo(t) + / 2 (07i,i(*.*) " +wu {t-r)+ &\t), 

bpi(t) - dp 2 (t) + e, 

2[a(f) 7 i,i(f, t) - c 7 i, 2 (t, t)} + 2w Ul (t - r) ( 1A {t, t - r) + (3 2 , 
2[67i,2(*,t) -d72, 2 (t,t)], 

bj 1A (t,t) + [a(t) - d]7i,2(M) " C72 )2 (M) +W* " T ) {*,i(t,t- r), 

(3 2 



2[a(t)p lil (t,t)-cp 1 , 2 (t,t)] + 2wu 1 (t-T)p 1 , 1 (t,t-T) + 



N ' 



2[bp li2 (t,t)-dp2, 2 (t,t)], 

bp hl (t,t) + [a(t) - d]p lj2 (t,t) - cp 2i2 (t,t) +wu 1 (t - r)p 2il (t,t- r), 



a(t) = f 1 (t) + 3f 3 (t) ll , 1 (t,t), 
u (t) = 9o(t) + 92(thi,i(t,t), 



(8) 
(9) 
(10) 

(11) 
(12) 

(13) 

(14) 

(15) 



(16) 
(17) 
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u 1 (t)=g 1 (t) + 3g 3 {t)j 1 , 1 (t,t), (18) 

CU*, t') = (jt=t[) \ n p^ t') - i.At, t% (19) 

where f t (t) = (l/e^F^^t)) and g t (t) = (l/£!)GW(pi W)- Equations (8)-(15) include the 
higher-order terms of j Ktl/ (t,t — r) and p Kjl/ (t,t — r), whose equations of motions are given 
by (m > 1) 

= [a(t) + a(t - mr)]7i 5 i(t, t - mr) - c[y 1:2 (t, t - mr) + 72,1 (i, t - mr)] 
+ w [u± (t — r) (t — r, t — mr) 

+ Ul {t-(m+ l)r) Ci,i(*» * - i m + i) T )} + P 2 AH, (20) 
= 6[7i,2(*, t - mr) + 72i i(t, i - mr)] - 2d 72i2 (t, t-mr), (21) 



d7i,i(*>* ~ 


mr) 


dt 






mr) 


dt 




071,2(4, * - 


mr) 


dt 




dl2 i(t, t — 


mr) 


dt 




dpi,i(t,t - 


mr) 


dt 




dp 2 ,2(t,t- 


mr) 


dt 




dp 1>2 (t,t - 


mr) 


dt 




dp 2 ,i(t,t- 


mr) 


dt 



= 671,1 t - mr) + [a(t) - d]^ li2 {t,t - mr) - C7 2)2 (t, t - mr) 

+ wu\ it - r) Ci, 2 it T,t mr) , (22) 

= bj lt i(t,t - mr) + [ait - mr) - d]^ 2 ,iit,t - mr) - 072,2 (M - mr) 

+ w Wl (t - (m + l)r) C 2 ,i(*, * - (m + l)r), (23) 

= [ait) + ait - mr)]pi i i(t, t - mr) - c[pi, 2 (i, t - mr) + p 2 ,i(t, t - mr)] 

+ w[ui(t — r)p 1)1 (t — r,t — mr) 

+ Ul (t-(^ m + l) T ) Pltl (t,t-(m + l)T)] + (^ A(mr), (24) 

= b\p lj2 (t, t - mr) + p 2 ,i(t, t - mr)] - 2dp 2j2 (t, t - mr), (25) 

= bpi tl {t,t - mr) + [ait) - d]p lj2 it,t - mr) - cp 2>2 it,t - mr) 

+ wui(t - r)p lj2 it - r,t - mr), (26) 

= 6pi i i(t, t — mr) + [a(t — mr) — d]p 2 ,i(i, t — mr) — cp 2j2 (t, t — mr) 

+ wui(t-(m+l)T)p2,i(t,t- (m+ l)r), (27) 



where A(rr) = 1 for x = and otherwise. 



C. Summary of our method 

The original two-dimensional SDDE given by Eqs. (1) and (2) are transformed to infinite- 
dimensional deterministic DDEs given by Eqs. (8)-(15) and (20)-(27), which are due to 
non-Markovian property of SDDE. It is, however, impossible to simultaneously solve these 
infinite-order recursive equations. We will adopt the level-m AMM (AMMm) in which the 
recursive DEs are terminated at the finite level m, as 
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IkA 1 -. t-(m + l)r) = 7 /S) „(i, * - mr), 
Pic,i/(<, * - ( m + = /»«,!/(*, * — (m+ l)r), 
5-i (t - (m + l)r) = 5-i (t - mr), 



(28) 
(29) 
(30) 



leading to 8(m + l)-dimensional DEs. In the following Sec. Ill, we will examine AMMm, 
performing calculations with changing m, in order to show that AMM5 may yield results 
in fairy good agreement with those of DS [Fig. 5(b)]. In the limit of r = 0, Eqs. (20)-(27) 
reduce to Eqs. (10)-(15), then Eqs. (8)-(15) agree with Eqs. (20)-(27) in Ref. [23] for FN 
neurons ensembles without delays [26]. 

Model calculations will be reported in the following Sec. III. DSs have been performed 
for 2N DEs given by Eqs. (1) and (2) by using the fourth-order Runge-Kutta method with 
a time step of 0.01. Initial values of variables at t G (— r, 0] are Xi(t) = Ui{t) = for % — 1 to 
N. DS results are the average of 100 trials otherwise noticed. AMM calculations have been 
performed for Eqs. (8)- (30) by using also the fourth-order Runge-Kutta method with a time 
step of 0.01. Initial values are fii(t) = /x 2 (0 = at t G [— r, 0], and 7 Kjl ,(t, t') = p K ^(t,t') = 
t G [— t, 0] or t' G [— r, 0] (t > t'). All calculated quantities are dimensionless. 



In this study, we pay our attention to the response of the FN neuron ensembles to a 
single spike input of I^ e \t) given by [23] 



where Q(x) = 1 for x > and otherwise, A stands for the magnitude, t{ n the input time 
and T w the spike width. We have adopted the same parameters of A = 0.10, t in = 100 and 
T w = 10 as in Ref. [23]. Parameter values of w, r, (3 and iV will be explained shortly. 

When an input spike given by Eq. (31) is applied, the oscillation may be triggered when 
model parameters are appropriate. The w-t phase diagram showing the oscillating (OSC) 
and non-oscillating (NOSC) states is depicted in Fig. 1, which is calculated for (3 = and 
iV = 10. In the case of f3 — 0.01, for example, the OSC region is slightly shrunk compared to 
that for (3 = 0, as will be shortly discussed [Figs. 5(a) and 5(b)]. The w-t phase is separated 
by two boundaries in positive- and negative--?/; regions. Circles in Fig. 1 express pairs of w 
and r adopted for calculations to be shown in Figs. 2 and 3. Along the horizontal, dashed 
line in Fig. 1, the w value is continuously changed in calculations to be shown in Figs. 4(a) 
and 4(b). 

In order to monitor the emergence of the oscillation, we calculate the quantity: 



III. MODEL CALCULATIONS 



A. Effects of coupling (w) and delay (r) 



I^ e \t) = A Q(t - t in ) 9(f in + T W — t), 



(31) 




(32) 



with 



o(t) 



-]T[< Xi (t) 2 > - < Xi (t) > 2 ], 

i 

M0 2 -^) 2 + 7i,i(0> 



(33) 



(34) 
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which becomes finite in the oscillation state but vanishes in the non-oscillating state, the 
overline denoting the temporal average between t\ (=2000) and £2 (=4000). 
The synchrony within ensembles is measured by [22] [23] 



°s = S(t), 



(35) 



with 



S(t) = 



N — 1 



(36) 



which is and 1 for completely asynchronous and synchronous states, respetively. 

We have calculated time courses of /^i(t), 71,1 (t, t), pi,i(t, t) and S(t), whose results are 
depicted in Figs. 2(a)-2(l), solid and dashed curves denoting results of AMM and DS, 
respectively. 

For r = 0, an output spike of pi(i) fires after an applied input which is plotted at the 
bottom of Fig. 2(a) [and also of 2(e) and 2(i)]. It is noted that state variables are randomized 
when an input spike is applied at t — 100 because independent noises have been added since 
t = 0. Figures 2(b) and 2(c) show 7^1 and p^i for r = 0, respectively. The synchronization 
ratio S(t) for r = shown in Fig. 2(d) has an appreciable magnitude: its maximum values 
calculated in AMM are 0.038 and 0.077 at t — 107 and 123, respectively. Figure 2(e) shows 
that when a delay of r = 20 is introduced, an input signal leads to a spike output with 
an additional, small peak in p,\ at t — 133. Figures 2(f) and 2(g) show that although a 
peak of 7i ; i for r = 20 becomes larger than that for r = 0, a peak of p\\ is decreased by 
an introduced delay. Maximum values of S(t) calculated by AMM are 0.154 and 0.130 at 
t = 126 and 140, respectively, for r = 20. We note from Fig. 2(i) that for a larger r = 60, 
an input spike triggers an autonomous oscillation with a period of about 65. Peaks in 7^, 
P11 and S are progressively increased with increasing t as shown in Figs. 2(j), 2(k) and 2(1): 
peaks of 7^1, p^i and S saturate at t ~ 1200 with the values of 0.00253, 0.00014 and 0.098, 
respectively, in AMM calculations. We note in Figs. 2(a)-2(l) that results of p\ obtained by 
AMM and DS are indistinguishable, and that AMM results of 7^, p ltl and S are in fairly 
good agreement with those of DSs. 

Figure 1 shows that although the obtained NOSC-OSC phase is nearly symmetric with 
respect to the w = axis, it is not in the strict sense. Actually the property of the 
oscillation for inhibitory couplings (w < 0) is different from that for excitatory couplings 
(w > 0). Figures 3(a) and 3(b) show autonomous oscillations for w — 0.1 and w = —0.1, 
respectively, with r = 60, (3 = 0.01 and N = 10. The period of the oscillation T is given 
by T = r + Tj where denotes the intrinsic delay for firings. For inhibitory feedback with 
negative w, FN neurons fire with the rebound process, which requires a larger Tj for firing 
than for excitatory feedback with positive w. Then the period of T = 86 for autonomous 
oscillation with the negative w becomes larger than that of T = 65 with the positive w. 

By changing the w value along the horizontal, dashed line in Fig. 1, we have calculated 
the w dependence of a and a s , whose results are plotted in Figs. 4(a) and 4(b), respectively, 
for (3 = 0.0001 and 0.01. The oscillation emerges for w ~ 0.058 or w ~ —0.063 with 
(3 = 0.0001, while with f3 = 0.01 it occurs for w ~ 0.060 or w ~ —0.070. The transition from 
NOSC to OSC states is of the first order because a Q is abruptly increased at the critical 
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coupling of w — w c , where a s has a narrow peak. In contrast, the relevant NOSC-OSC 
transition in the nonlinear Langevin model is of the second order [22]. 

We have investigated, in more detail, the w dependence of a a and a s near the transition 
region of 0.05 < w < 0.07, which is sandwiched by vertical, dashed lines in Figs. 4(a) and 
4(b), results for (3 = 0.0001 and (3 = 0.01 being plotted in Figs. 5(a) and 5(b), respectively. 
Figure 5(a) shows that the critical w value for the NOSC-OSC transition is w c ~ 0.0579 for 
(3 = 0.0001 both in DS and AMM5. When we adopt AMM1, we get the result showing the 
NOSC-OSC transition at w ~ 0.6, although we cannot get solutions for 0.0586 < w < 0.060. 
With the use of AMM2, we get the transition at w ~ 0.058, though solutions are not 
obtainable for 0.0580 < w < 0.0582. We have noted that AMMm converges at the level 
m — 3, above which calculated results are almost identical. Figure 5(b) shows that the 
critical value of w c for (3 = 0.01 is 0.0600 in DS and 0.0607 in AMM5. For m = 1, 2 and 
3, the NOSC-OSC transition occurs at w = 0.0644, 0.0609 and 0.0807, respectively: w c for 
m = 3 approaches that for m = 5 (in what follows results of AMM5 will be reported). 
It is interesting to note in Figs. 5(a) and 5(b) that the synchrony a s shows fluctuation- 
induced enhancement at the NOSC-OSC transition. This is due to an increase in the ratio 
of pi,i(£, t)/7i i i(t, t) in Eq. (36) although both pi i i(t,t) and 7i,i(t,t) are increased at the 
NOSC-OSC transition. Similar phenomenon has been reported in the nonlinear Langevin 
model [22] and in heterogeneous systems in which the oscillation emerges when the degree 
of the heterogeneity exceeds the critical value [33] [34] . 



B. Effects of noise (/?) 

Comparing Fig. 5(b) with Fig. 5(a), we note that when the noise intensity is increased 
form (3 = 0.0001 to (3 = 0.01, the critical w c value for the NOSC-OSC transition is increased: 
w c =0.0579 (0.0579) for (3 = 0.001 and w c =0.0600 (0.0607) for (3 = 0.01 in DS (AMM). Figure 
6(a) shows the (3 dependence of a Q and a s for r = 60, w = 0.06 and N = 10. a Q is rapidly 
decreased at (3 ~ (3 C where a s has a broad peak: (3 C is about 0.01 in DS while it is about 
0.0075 in AMM. Figure 6(b) shows that the similar (3 dependence of a and a s is obtained 
also for a larger w = 0.062, for which f3 c ~ 0.015 in DS and (3 C ~ 0.014 in AMM. A 
suppression of the oscillation by noises is realized in the Langevin model [22] and in some 
calculations for systems with heterogeneity [34], although the noise-induced oscillation is 
reported in Refs. [21] [35] [36]. In particular, Zorzano and Vazquez [21] (ZV) showed the 
noise-induced oscillation in FN neuron ensembles (N = oo) with time delays by using FPE 
method. The difference between ZV's results and ours may be due to the difference in the 
adopted FN model and/or ensemble size. In order to get some insight on this issue, we have 
performed AMM calculations for our FN model with larger ensemble sizes of N = 100 and 
1000, and obtained again a suppression of the oscillation by noises [37]. It is not clear for us 
how ZV took into account the non-Markovian property of SDDE within their FPE method 
[3] [5]. 
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C. Effects of size (TV) 



The N dependence of a Q and a s for (3 = 0.01, w = 0.06 and r = 60 is shown in 
Fig. 7 where open circles (squares) express a Q (a s ) in DS, and where thin (bold) solid curves 
denote a Q (a s ) in AMM. It is shown that with increasing the size of ensemble, a Q is gradually 
increased at N ~ iV c where a s has a broad peak, the critical dimension being N c ~10 in DS 
and N c ~ 100 in AMM. Results of our AMM calculations are qualitatively similar to those 
of DS although calculated N c values are different between the two methods. 

IV. CONCLUSIONS AND DISCUSSIONS 

In Sec. II, we have obtained the infinite-dimensional ordinary differential equations. It 
is, however, possible to get expressions given by partial differential equations (PDEs) if we 
define the correlation functions: 

c KtX (t, z) = ^J2< Sx M 6x ^( f - z ) >> ( 37 ) 

i 

D K>x (t, z) = < 5X K {t) 5X x (t - z) >, (38) 

introducing a new variable z [see Eqs. (5) and (6)]. For example, PDEs for Ci^{t,z) are 
given by 

9Chl d f 0) = 2[aC 1>1 (f, 0) - cC 1;2 (t, 0)] + 2w Ul (t - r)E 1A (t, r) + [3 2 , (39) 

+ wu 1 (t-T)E 1:1 (t-T,z-r), for^>0 (40) 

where E hl (t,z) = [ND 1A (t,z) - C 1A {t, z)]/[N - 1]. It is noted that Eqs. (39) and (40) 
correspond to Eqs. (10) and (20), respectively. Then we have to solve PDEs including 
fJ> K (t), C K! x(t,z) and D K \(t,z) with a proper boundary condition in the (t,z) space. A 
similar PDE approach has been adopted in Ref. [6] for an analysis of the stationary solution 
of the linear Langevin equation with delays. In an earlier stage of this study, we pursued 
the PDE approach. We realized, however, from the point of computer programming that 
the use of the ordinary DEs given in AMM is more tractable than that of PDEs. 

Our calculations have shown that FN neuron ensembles with delays exhibit the multi- 
stability when model parameters such as w, r, f3 and iV are varied. The multistability is the 
common property of the system with time delay. Actually the nonlinear Langevin ensem- 
bles discussed in I also show the multistability: the w — r phase diagram of FN ensembles 
shown in Fig. 1 is similar to that of the Langevin ensembles shown in Fig. 6 of I. In either 
case, fluctuation-induced synchronization is realized near the transition between OSC and 
NOSC states. These results imply that the oscillating, highly synchronous state may be 
realized in ensembles for smaller couplings with a proper delay than with no delays. This is 
consistent with the recent result of Ref. [38] , where the importance of delays is stressed for 
the long-range synchronization with low coupling strength. 
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In summary, we have discussed dynamics of FN neuron ensembles with delays by using 
a semi-analytical method developed in I. Our method has a limitation of weak noises but 
it is free from the magnitude of delay times. This is complementary to the small-delay 
approximation [3], whose application to FN neuron ensembles with delays is discussed in 
appendix C. For FN ensembles to show the oscillation, we have to adopt an appreciable 
magnitude of delay (r ~ 20), for which SDA method cannot be employed. In this study 
we have discussed only the case of a single spike input. Our method may be, however, 
applicable to arbitrary inputs such as periodic spike trains and Poisson spikes, as was made 
for HH neuron ensembles (without delays) [24] . Although results calculated by our method 
are in fairly good agreement with those obtained by DC, the quantitative analytical theory 
is still lacking. In this study, we have assumed regular couplings (wij = w) and uniform time 
delays (r^ = r). In real systems, however, couplings are neither regular nor random, and 
time delays are nonuniform with a variety of dendrite radius and length. It is interesting to 
include these properties by extending our approach, which is in progress and will be reported 
in a future paper. 
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APPENDIX A: DERIVATION OF EQS. (8)-(15) 

We express Eqs. (1) and (2) in a Taylor expansion of 5xi (= Sxu) and Syi (= 5x 2 i) up to 
the third-order terms to get 

= h(t)5x t (t) + f 2 [5xi(t) 2 - 7i,i(M)l + f3(t)6 Xi (t) 3 - c5y t (t) 
+ Ut) + Slt\t-r), (Al) 
dSyi ^ =b6xi(t)-d6 yi (t), (A2) 



dt 
with 



5l\ c \i) = w 



£ + #T E N«f - 7ul + #T E • (A3) 



where f e (t) = (l/£ ^^(^(t)) and g e (t) = {1 / £ \)G^ (^(t)) . Averages of Eqs. (Al) and 
(A2) with Eqs. (3) and (4) yield DEs for means of d\i\jdt and d\x 2 jdt [Eq. (8)]. DEs for 
variances and covariances may be obtained by using the equations of motions of 5xi and Syi. 
For example, DE for d^i :2 (t, t)/dt is given by 

^4?<(^W+M^h (A4) 

which leads to Eq. (12). DEs for other variances and covariances are similarly obtained. 
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APPENDIX B: DERIVATION OF EQS. (20) AND (27) 

In the process of calculations of Eqs. (8)-(15), we get new correlation functions given by 

Si(*i,* 2 ) = ^E < Uh) >, (Bi) 

i 

S 2 (h,t 2 ) = ^E < <%(*i) €i(* 2 ) >, (B2) 

i 

where <5xj = Sxu, 5yi = 5x 2 i, h — t and t 2 — t — mr, or ti — t — mr and t 2 = t. We will 
evaluate them by using DEs for Sxi(t) and Syi(t), which are linearized from Eqs. (A1)-(A3): 

^£ = a(t)5x t (t) - cS Vi (t) + (j^) E 9i(t~ r)5x 3 {t - r) + (B3) 

^ML = bSx i (t)-dSy i (t) 1 (B4) 

where a(i) = /i(i) + 3/ 3 (t)7i ) i(t, t). Neglecting the t dependence in a(t), we get formal 
solutions of Eqs. (B3) and (B4) given by 

A + d\ ft 



6xi(t) = [j^g) J ds ^ {t ' S)A i(w=l) E9i(s - r)5x 3 (s - r) + Us)] 

dsexp (ts)B { ,__} y gi{s _ T)6X]{s _ T)+Us)]} m 



Syi{t) 



A-B 



A 



w \ 


N — 


ij 


w 




N — 


t) 


w 


i) 


N — 




f w 


t) 


\N- 





Mi) 



Mi) 



(B6) 



where A and B are roots of the equation given by z 2 — (a — d) z — a d + b c = 0. By using 
the method of steps in Ref. [6], we obtain the step by step functions, from which we get 

S l {t,t-mT) = S 1 {t-mr,t)= ^y^J A(mr), (B7) 
S 2 {t,t-mr) = S 2 {t-mr,t) = 0, (B8) 

where A(x) = 1 for x = and otherwise. By using Eqs. (B7) and (B8), we get Eqs. 
(20)-(27). The assumption of a neglect of the t dependence in a(t) may be justified, to some 
extent, from results calculated by our method which are in fairly good agreement with those 
by DS as reported in Sec. III. 



APPENDIX C: THE SMALL-DELAY APPROXIMATION 

When the delay r is very small, we may adopt the small-delay approximation (SDA) 
proposed in Ref. [3]. With this approximation, we first transform the SDDEs to stochastic 
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non-delayed DEs, and then to deterministic DEs with the use of DMA [23]. For a small r, 
we may expand Xu(t — r) in Eq. (1) as 



xu(t -r) ~ xu(t) - r 



dxuit) 
dt 



(CI) 



with which Eq. (1) becomes stochastic non-delayed DEs given by 



dx u (t) 
dt 



+ (j^j) E G\x l3 {t)) 



dxij (t) 
dt 



= F{x u ) - cx l2 + (-^—) £ G( Xlj (t)) + m + I (e) - 



(C2) 



When we apply DMA to 2iV-dimensional stochastic DEs given by Eqs. (2) and (C2), we get 
equations of motions for means, variances and covariances, given by 



= [1 - wr Ul ][f (t) + / 2 (07i,i(*, t) - cp 2 (t) + wg (t) + /<«>(*)], 



dp 2 (t) 
dt 

ri7i,i(M) 
dt 



bfii(t) - rf/i 2 (t) + e, 

2[a(*)7i,i(^ " C7i, 2 (t, t) + w Ul (t)Ci,i(t, t)} + 1 



— 2wrui(t) 



WUi(t) 



^72,2 (M) 

dt 

dji t2 (t,t) 
dt 



dpi,i{t,t) 
dt 

df)2,2{t,t) 

dt 

dp li2 (t,t) 
dt 



a(t){ ltl (t,t)-ct li2 (t,t) + [j^y J (Npi,i(t,t) -Ci,i(*>*)) 
2[&7i,2(t,t)-d72,2(M)]> 

&7i,i(*, *) + [«(*) " ^]7i, 2 (t, - C7 2 , 2 (t, t) + umi(*)Ci,2(*, *) 



-WTUi(t) 



a(*)Ci,2(*, - <2, 2 (^ + (^Tj) (Wpi )2 (*, - Ci, 2 (t, t)) 



2[1 - wrui(t)] 



a(t)pi,i(i, i) - cpi j2 (t, i) + wwi(t)pi i i(t, t) + 



2N 



2[b Pl , 2 (t,t) - dp 2t2 (t,t)}, 

bp hl (t,t) + [a(t) - d]p lj2 (t,t) - cp 2t2 (t,t) + wui(t)pi )2 (t,t) 
-WTUi(t)[a{t)p 1)2 (t,t) -cp 2>2 (t,t) +wu 1 (t)pi j2 (t,t)], 



(C3) 
(C4) 

(C5) 
(C6) 

(C7) 
(C8) 
(C9) 

(CIO) 



where a(£) and ( Kt \(t,t) are given by Eqs. (16) and (19), respectively. 

A numerical comparison between AMM and SDA is made in Fig. 8, where solid and 
chain curves denote results of AMM and SDA, respectively. For r = both methods lead 
to the identical result. For small delays of r = 1 and 2, results of SDA are in fairly good 
agreement with those of AMM. As the delay is increased to r > 5, however, the discrepancy 
between the two methods becomes significant. 
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FIGURES 



FIG. 1. The w-t phase diagram showing the oscillating (OSC) and non-oscillating (NOSC) 
states for (3 = and N = 10. For sets of parameters of w and r marked by circles, time courses of 
pi(t), j(t,t), p{t,t) and S(t) are calculated, whose results are shown in Figs. 2 and 3. Along the 
horizontal dashed line (r = 60), the w dependence of a Q and a s is calculated in Figs. 4 and 5. 

FIG. 2. (color online). Time courses of p>i(t), 71,1 (i), Pi,i(t) and S(t) calculated by AMM 
theory (solid curves) and DS (dashed curves) with A = 0.10, (3 = 0.01, w = 0.1 and N = 10: (a) 
Ml, ( b ) 7i,i, ( c ) Pi,i and ( d ) 5 for r = 0, (e) fii, (f) 71,1, (g) p ltl and (h) S for r = 20, and (i) m, 
(j) 71,1, (k) pi t i and (I) S for r = 60. Chain curves at bottoms of (a), (e) and (i) express input 
spikes. 

FIG. 3. Time courses of p\{t) showing the oscillation for (a) w = 0.1 and (b) w = —0.1 with 
r = 60, (3 = 0.01 and N = 10 calculated by AMM, the result of (a) being shifted upwards by 2. 

FIG. 4. The w dependence of (a) a Q and (b) a s for (3 = 0.0001 (solid curves) and (3 = 0.01 
(dashed curves) with r = 60 and iV = 10. The region sandwiched by dashed, vertical lines is 
enlarged in Figs. 5(a) and 5(b) for (3 = 0.0001 and 0.01, respectively. 

FIG. 5. The w dependence of a a and a s for (a) (3 = 0.0001 and (b) (3 = 0.01 with r = 60 and 
N = 10. Thin and bold solid curves denote results of 10 a Q and a s , respectively, in AMM, whereas 
squares and circles express those of 10 a and a s , respectively, in DS. AMM results with different 
level m (=1, 2, 3 and 5) are shown. Dotted lines are only for a guide of the eye (see text). 

FIG. 6. The (3 dependence of a Q and a s for (a) w = 0.60 and (b) w = 0.62 with r = 60 and 
N = 10. Thin and bold solid curves denote results of 10 a a and a s , respectively, in AMM whereas 
squares and circles express those of 10 a Q and a S} respectively, in DS. Dotted lines are only for a 
guide of the eye. 

FIG. 7. The iV dependence of a Q and a s for (3 = 0.01, r = 60 and w = 0.06. Thin and bold solid 
curves denote results of 10 a and a s , respectively, in AMM, whereas squares and circles express 
those of 10 a Q and a s , respectively, in DS. Dotted lines are only for a guide of the eye. 

FIG. 8. The time course of p.i(t) calculated in AMM (solid curves) and in a small-delay ap- 
proximation (SDA) (chain curves) with (3 = 0.01, w = 0.1 and N = 10 (see appendix C). 
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